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We analyze an optical parametric oscillator (OPO) in which cascaded down-conversion occurs inside a cavity 
resonant for all modes but the initial pump. Due to the resonant cascade design, the OPO present two x level 
oscillation thresholds that are therefore remarkably lower than for a x' 3 ' OPO. This is promising for reaching the 
regime of an effective third-order nonlinearity well above both thresholds. Such a \^ cascaded device also has 
potential applications in frequency conversion to far infra-red regimes. But, most importantly, it can generate 
novel multi-partite quantum correlations in the output radiation, which represent a step beyond squeezed or 
entangled light. The output can be highly non-Gaussian, and therefore not describable by any semi-classical 
model. 
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I. INTRODUCTION 

Continuous-variable (CV) quantum information is an in- 
teresting flavor of quantum information (QI) fljj, |2|]. While 
easily implemented by use of well established quantum op- 
tics techniques, benefiting from large flow rates and broad 
spectral band, it has long been based on coherent states and 
linear Bogoliubov transformations (quadratic Hamiltonians) 
and therefore restricted to positive Wigner functions of Gaus- 
sian character. These states are not general enough for uni- 
versal quantum information operations [3]. For instance, it 
has been shown that quantum computation based solely on 
Gaussian CV states can be efficiently simulated by a classical 
computer 0]. Also, CV entanglement purification requires a 
Kerr-nonlineari ty -based QND measurement |5] or, in general, 
a non-Gaussian state @]. However, it has also been shown 
that one-way quantum computing can be implemented using 
Gaussian cluster-state entanglement combined with non-CV 
(e.g., photon counting) measurements Q]. 

Recently, successful "degaussification" experiments, using 
homodyne detection conditioned on single-photon detection 
IE |H [T(J, [Till , have successfully generated negative Wigner 
functions from initial squeezed states. Here, we investigate 
different type of sources, which can produce non-Gaussian 
light directly. Theoretical studies of optical parametric oscil- 
lators (OPO), which are based on a single second-order optical 
nonlinearity (y^) have shown non-Gaussian signatures to be 
rather scarce 111211 except in the case of the tripartite correla- 
tions between the three fields [13]. An interesting approach 
is to use an optical nonlinearity of, at least, third order. This 
has been theoretically investigated ITil [l5l [T6L \vh . In prac- 
tice, a based OPO would have the problems of requiring a 
very large and possibly prohibitive input power threshold for 
downconversion, together with an even higher threshold for 



the onset of nonclassical effects, such as the formation of star 
states 0]. 

In this paper, we show how the use of a cavity- resonant cas- 
cade of second-order nonlinearities can yield a low-threshold 
OPO which possesses the effective behavior of a \^ OPO in 
certain regimes and is more accessible experimentally. Note 
that related systems have been studied before, in the purely 
classical case and for completely different purposes, such as 
producing new tunable optical sources in the infrared 111 811 or 
achieving optical phase-locking in a 3:1 frequency ratio for 
frequency metrology |[T9ll20ll . Parametric amplifiers and os- 
cillators have become a widely used, even standard part of the 
repertoire of laser physics and quantum optics MM- Above 
the classical threshold points, these devices are a useful tool 
for frequency conversion. Below threshold, quantum effects 
dominate, leading to squeezing and entanglement. These de- 
vices that rely on non-resonant, nonlinear optics interactions 
have proved experimentally superior to other resonant or near- 
resonant alternatives, due to the fact that absorption is sup- 
pressed. 

There are other possible quantum effects available, as well 
as direct down-conversion in the linear regime well below 
threshold. For example, exploration of non-equilibrium quan- 
tum criticality is possible near threshold. This results in large 
critical fluctuations and phase-transitions. The fluctuations in 
this case become non-Gaussian, but the dominant critical fluc- 
tuations have a rather classical character. Here, we explore an- 
other path to such non-Gaussian behavior, in which extremely 
nonclassical correlations are generated through the presence 
of a second down-conversion crystal placed inside the cavity. 
We show that this results in an intricate pattern of new phase- 
transitions at the classical level, in which there are two distinct 
threshold points. At the quantum level, below the first thresh- 
old, there are very strong triple correlations between the three 
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down-converted modes, which have no classical analog. 

This paper is structured as follows. In Section II, we ex- 
plain the basic model and the theoretical phase-space tech- 
niques that are used here. In Section EH, we present an analyt- 
ical study of the system's stationary solutions. In Section IV, 
we turn to a treatment of stability properties and fluctuations 
in one particular type of down-conversion scenario. In Section 
V, we give numerical simulations of more general cases which 
can also yield regimes of interest. These simulations demon- 
strate the stability regions, in the same spirit as was achieved 
for the x {2) OPO JH. We give conclusions in Section VI. 




Figure 1: (Color online) Schematic of the resonant cascaded down- 
conversion system. A driving field, that is pumped at a frequency cjq 
with amplitude Eo, enters a cavity that contains two x nonlinear 
crystals. The first xi crystal down-converts the original mode ao 
into two modes oi and a 2 with frequencies ui\ and W2, respectively. 
Then the second mode &2 undergoes a further down-conversion, via 

(2) 

the \ 2 crystal, into the two modes £3 and 0,4 with frequencies ^3 
and W4, respectively. 



II. ANALYTICAL TREATMENT 
A. Cascaded parametric oscillator model 

The model system for cascaded down-conversion consists 
of two quadratically nonlinear elements with nonlinearities xi 
and \2 inside an optical cavity (c.f. Fig.Q]). The cavity sup- 
ports five resonant modes at frequencies u>i (i = 0,1, ...4). 
The mode loq is the pump mode, driven by an external coher- 
ent driving field at the same frequency ujq. The cavity modes 
u>i are described by creation and annihilation operators a] and 
dj with commutation relations [Si,aj] = 5ij. The first non- 
linear element converts the pump mode uio into the signal and 
idler modes ui\ and L02 by means of nondegenerate parametric 
down-conversion, where uiq = lo\ + 0J2 (u>i 7^ lo-i)- The sec- 
ond nonlinear crystal supports down-conversion of the mode 
u>2 into the second pair of signal and idler modes, U3 and 0J4, 
where u>2 — ^3 + lo^. We will call the field S2 at u>2 the "in- 
termediate pump." The modes may decay via cavity losses at 
the respective rates 7,;, i € [0, 4]. 

In the absence of the optical cavity, this interaction consti- 
tutes a cascade of quantum systems in the sense investigated 
by several authors before [23], where the second stage does 



not feed back to the first stage. Here, the situation is differ- 
ent precisely because of the cavity feedback, hence our use 
of the term resonant cascade throughout the paper. Within 
this frame, we will distinguish two situations: the first one is 
the nondegenerate resonant cascade, for which the fields Si, 
S3, and CI4 are distinguishable (i.e., lo\ 7^ UJ3 7^ u>4, or hav- 
ing different polarizations or wave-vector directions). In this 
case, the only physical observable affected by both stages of 
the cascade is the intermediate pump a%. This is the case that 
will be investigated analytically, with additional simplifying 
hypotheses, and numerically, without those hypotheses. The 
second case is the degenerate resonant cascade, for which the 
signal fields are indistinguishable: Si = S3 = S4 (and hence 
u>i = 0J3 = L04), In that case the signal field and the interme- 
diate pump interact in both nonlinear media and the dynamics 
are richer. That case will be explored by numerical simula- 
tions. Obviously, intermediate situations do also exist, e.g., 
u>i = 0J3 7^ 0J4, but we will not consider them here. 

B. Hamiltonian and equations of motion 

The model Heisenberg-picture Hamiltonian for the system, 
in the rotating-wave approximation, is given by: 

4 

H = S^iStSj + ih(E e- luJot al - E*e iuJot a ) 

i=0 

+ ihxi(a a\al — SjSiS 2 ) + zfr^^SgSj — S2S3S4) 
4 

+ ^(S i fl+Stf i ) (2.1) 

i=0 

Here, Eq describes the complex amplitude of the driving field. 
The coupling constants xi an d X2 are proportional to the 
second-order susceptibilities of the two nonlinear elements, 
respectively. We assume that they are positive, without loss 
of generality, since phase factors can always be absorbed into 
the definitions of the mode functions and their operators. The 
operators f , and f J describe the coupling of each intracavity 
mode to the reservoir of external modes. These give rise to 
the losses of the cavity modes u>i at rates 7.;. 

1. Master Equation 

Transforming to an interaction picture in which all opera- 
tors are transformed to rotating frames, i.e., 

aj(t) = Sj-e" 1 ^ 4 , (2.2) 

one can derive the following master equation for the system 
density operator p i24ll : 

^ = [Eoal - E^a ,p] + Xi[SoS|sJ - SjSiS 2 ,/5] 

+ X2[a2S^S 4 - S^S 3 a 4 , p] 
4 

+ ^2 7i(2SipSj - pa\a % - aJSj/3). (2.3) 

i=0 
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While in principle this master equation can be solved numer- 
ically in a number-state representation, in practice this is not 
possible. The complexity of the Hilbert space — especially 
for this five mode problem — is enormous, given any moder- 
ate number of photons present in the five interacting modes. 
Instead, we solve this problem using phase-space representa- 
tion methods, such as the positive-P representation[25|. 

C. Positive-P representation 



together with the corresponding equations for . Here, the 
dots imply a time derivative, and the terms Q (t) and £ 4 - (t) are 
independent complex Gaussian noise sources with zero means 
and the following nonzero correlations: 



(Ci(W)) = <d + (i)C 2 + (i')) = S(t t'), 

(muv)) = (ct(t)ct(t')) = s(t - 1). (2.6) 



Using the positive-P representation we can transform the 
master equation, Eq. (12.3b . into a Fokker-Planck equation 02511 
expressed as: 

d 

-(7oao - E + xiaia 2 ) 







dt 



da 

d . +x 

+ « — (71 "l ~ Xi^oaJ) 
a oi\ 

9 , 4- x 

+ -Z — (72«2 - Xiaoai + X2«3a4j 
oa 2 

+ q (73"3 - X2a2«4 ) 

oa 3 

+ T. (74 "4 - X2a2«3 ) 

OCX4 



C) 2 



+ 



da\doi2 

d 2 
da 3 da^ 



(Xiao) 

(X2tt2) + h.c. 



P(a,a + ,t). 

(2.4) 



Here, a = 



(ao, ai, «2, CK3, 014) and a + = 
) represent the sets of coherent 
state amplitudes a, and af in the expansion of the density 
operator in terms of the positive P-representation, corre- 
sponding to the annihilation and creation operators dj and a\ . 
We recall that in the positive P-representation, the amplitudes 
on and af are independent complex c-numbers, and h.c. 
in Eq. ( 12.41) represents the terms equivalent to Hermitian 
conjugate operators, obtained from the previous terms by 
replacing aj — » a+ and vice versa, while Eq is replaced by 
Eq . The transformation requires an assumption of vanishing 
boundary terms which can be checked numerically. This is 
generally extremely well-satisfied ll26Tl for these open systems 
provided Xi << 7j> which is typically the case in nonlinear 
optics experiments. If required, further stochastic gauge 
transformations [27] can be used to eliminate boundary terms. 

The Fokker-Planck equation ( 12.41 ) is equivalent to the fol- 
lowing set of stochastic differential equations I2H1 . in the Ito 
form: 

= -7o"o + E Q - xiaia 2 , 
a\ = -71 ai + Xia Q a^ + VXi"oCi(*), 
«2 = -72«2 + Xiao^t - X2"3a4 + VXi a oC2(*)> 
a 3 = -730:3 + X2a 2 aj + ^/x2^C 3 (t), 
ola — 74^4 + X2a 2 a 3 + N /X2"2C4(*)- (2.5) 



The above set of the stochastic equations of motion, Eq. d2.5l) , 
can be solved either numerically or else using approximate an- 
alytic treatments such as perturbation expansions around sta- 
ble semi-classical steady states. Quantum mechanical observ- 
ables that are expressed in terms of normally ordered opera- 
tor moments ^(d|) n (di) m ^) correspond to stochastic averages 
((a,)" 1 (ajr). 



D. The semi-classical theory 



We can also transcribe the master equation, Eq. ( 12.31 ). as 
a c-number phase space evolution equation using the Wigner 
representation 12911 



Pyy(a, a* 



d zxw(z,z*)e~ 



(2.7) 

where xs( z , z *)> the characteristic function for the Wigner 
representation , is given by 



Xw (z,z*)=Tr [pe 1 



(2.8) 



This transcription is particularly useful for semi-classical 
treatments in which we include quantum noise terms from 
the reservoirs, but neglect higher-order quantum noise from 
the nonlinear couplings. This approximation is also called a 
truncated Wigner approximation, as it is obtained from a full 
Wigner-Moyal equations via truncation of third-order deriva- 
tives. 

The equation for the Wigner function for the nondegenerate 
parametric amplifier that corresponds to the master equation 
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given by Eq. ( 12.31 ) turns out to be 1I23I1 



at 



(70^0 + Xiai«2 - £) 



\ — 

[da 

+ da* <a ° a *° + Xia * lOL * 2 ~ ^ 



d_ 

dai 
d 



d 

t^- (7i«i - Xia* 2 a ) + — (7ia x - XiO 2 a ) 



+ t; — (72^2 - Xi a i a o + X2O3O4) 
aa 2 

d 

+ ^— - (720-2 - XiOi«5 + X203Q;4) 

a (7303 - X2 4 a2j + (73O3 - X2040 2 



da 

d 



d 



+ (7404 - X2O3CX2) + (74O4 - X2O3O2) 



da 
+70 

+73 



(9 2 



d 2 



c) 2 



da da*, 

d 2 



da^da 



■71 
•74 



daicta^ da 2 da\ 
d 2 



daedal 

Q3 Q3 

+ 



Xi 

4 \da1da2daQ da\da*,dav 
X2 ( d 3 d 3 



+ 



4 V da^da^dat datda\da.2 



Pw 



It is common to drop the third order derivative terms, in an 
approximation valid in the limit of large photon number. This 
allows one to equate the resulting truncated, positive-definite 
Fokker-Planck equation with a set of stochastic equations. 



These are: 






do = 


-7000 - 


h E - xiOio-2 + y/7oVo(t), 


«i = 


-7XQ!! - 


I-Xioooa +V7i r ?i(*) ! 


OL% = 


— 7202 - 


h Xl«0Oi - X203«4 + \Zj2lT2(t), 


03 = 


-7303 - 


h X20 2 04 + ^/j3~m(t), 


d 4 = 


-7404 - 


I-X2O2O3 + y/j4T)4,(t). 



(2.9) 



together with the corresponding equations for af . Here, the 
conjugate equations have conjugate noises as in a normal clas- 
sical phase-space. The terms r)i(t) are complex Gaussian 
noise sources with zero means and the following nonzero cor- 
relations: 

(m(t)v*(t')) = s t] s(t-t>), 

If we compare the two sets of Ito stochastic equations, we 
see that the noise terms in the positive-P equations, Eq. ( 12.51 ). 
depend on the nonlinear coupling constant, while those in the 
Wigner representation, Eq. (12.9b . do not. 

The truncated Wigner theory can be regarded as a kind of 
hidden-variable theory, since it behaves as though the non- 
commuting quadrature variables were simple classical ob- 
jects. These equations imply that (aia\) = (hi) = 1/2 when 
there is no driving and no coupling, which is an expected 



result in a symmetrically-ordered representation. However, 
the truncation neglects third-order derivative terms which are 
present in the full Wigner equation, and are not always neg- 
ligible. The full Wigner theory is equivalent to quantum 
mechanics, and has no such limitations but it is no longer 
positive-definite, and therefore has no equivalent stochastic 
formulation. The advantage of the positive-P method is that it 
is able to generate stochastic equations without requiring this 
questionable truncation approximation. 



III. CLASSICAL STEADY STATES 

We first analyze the classical steady states of the system 
and then give the results of the linearized fluctuation analysis 
for their stability in the next section. In the classical limit, all 
quantum noise terms are neglected. The positive-P stochas- 
tic variables on and af are replaced by deterministic ampli- 
tudes a,; and a*, where a* is the complex conjugate of a^ 
and Eq. (12.5b then becomes 

do = -7000 + E — xiOia 2 , 

01 = — 7iai + xia a2, 

02 = -72O2 + XlOooJ - X20 3 «4, 

03 = -7303 + X20204, 

d 4 = -7404 + X20203- (3.1) 

The steady state solutions a° are obtained from Eqs. ( 13.11 ) by 
putting all time derivatives equal to zero, i.e., on — 0. We only 
consider the steady-state solutions in which (a ( -) + = (a®)*, 
as these correspond to classical fields. This corresponds to 
neglecting the effects of quantum fluctuations and consider- 
ing the equations for the mean field amplitudes a; = (a;), 
assuming that higher-order correlations factorize. 

The stability of the classical steady states with respect to 
small fluctuations can be checked by deriving the linearized 
equations of motion for the fluctuations Sai(t) = ai(t) — a° 
and Saf(t) = af(t) — (af)*. The steady states are stable 
provided all the eigenvalues of the appropriate drift matrix of 
the linearized equations have negative real parts. Here, we 
assume the following matrix form of the deterministic part of 
the linearized equations of motion: 



= Ax 



(3.2) 



where A is the drift matrix, and x denotes a column vector 
for fluctuations {<5a.;, Saf}. If the linearized eigenvalue anal- 
ysis reveals eigenvalues with non-negative real parts, this im- 
plies that the steady states are unstable. In this case, the lin- 
earized treatment of fluctuations around the classical steady 
states cannot be employed, and the equations of motion have 
to be treated exactly. 

To simplify our analysis and make analytic solutions avail- 
able, we will assume that the damping rates for all modes ex- 
cept the pump mode are equal to each other, 



7i = 72 = 73 = 74 = 7, 



(3.3) 
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Figure 2: (Color online) Steady state solutions for the scaled inten- 
sities rii/no_ cr of different modes (i = 0,1,. ..,4) as a function 

\E f /\E thrA \ 2 , for 



of the driving field intensity parameter e 2 

7o/t = 10. Here, Tio iCT = "f 2 /x 2 is me critical value of n[] at the 

Cthr-.i = l)- The second thresh- 
(I + 7/70) 2 = 1.21. 



first threshold, i?o| 2 = \Ethr,i\ 2 
old here corresponds to e 2 hr 2 



steady states (see Sec. IV) reveals that they are stable for driv- 
ing field intensities below a certain critical (threshold) value, 



\E f < \E thrtl \\ 



where 



\E 



thr,l I 



2 2 

,2 _ 7o7 
X 2 



(3.8) 



(3.9) 



is the first threshold. This allows us to introduce a dimension- 
less relative driving field parameter, 

_ \Eo\ 



\E t hr,l\ 

Thus, the first regime corresponds to conditions where both 
nonlinear crystals operate in the below threshold regime of 
parametric down-conversion. Here the steady state solutions 
for the modes loq, u>\ and u>2 are the same as in the usual non- 
degenerate parametric down-conversion with a single crystal 
IBfJl . Figure [2] plots the steady state solution jiq in the be- 
low threshold regime where we have also introduced a new 

|2/„,2 _ 7 2/ % 2^ Th j s j sme 



variable, namely n cr = | E thr 1 1 2 /i^ 



critical value of rig at the first threshold. 



while the pump mode is strongly damped, 



70 >> 7, 



(3.4) 



to model an interferometer which is not resonant at the pump 
wavelength. For simplicity we suppose that the coupling con- 
stants xi an d X2 are also equal: 



Xi = X2 = X- 



(3.5) 



Analysis of the equations of motion for the classical steady 
states reveals three different types of solutions, corresponding 
to three regimes of operation. 



2. First above-threshold regime 

In the first above-threshold regime, the amplitudes of the 
modes W3 and U14 remain zero, while the amplitudes of the 
pump, signal and idler modes (uiq, u>i and ui 2 ) are nonzero. 
Accordingly, we again use the intensity and phase variables, 
? exp(z<^) for i = 0,1,2, and write the 



0, 



and 



AO „,0 



steady state solutions as: 



1. Below threshold regime 

Here, the amplitudes of all intracavity modes except the 
pump mode ujq are zero, and we find that 



a\ = a° 2 = a% = a\ = 0, 
a° - ^ 

70 



(3.6) 



The last equation can be rewritten in terms of the steady state 
intensity n® = I o§ I (in photon number units) and phase 4>q 
(where = \fr^ cxp(i0q)): 

\E Q \ 2 hl 



4 



(3.7) 



where ipo is the phase of the driving field, i.e., Eq = 
\Eq\ exp(it^o). The linearized stability analysis of these 



n? 



o_ \E \ 

2 — 

X 

4>a = <^o, 



7o7 

i/2 ' 



(3.10) 



(3.11) 



We see that the steady state intensities n\ and nP, correspond 
to physical solutions (rt° > ) if the driving field intensity 
is above the first threshold, |i?o| 2 > \Ethr.i\ 2 - On the other 
hand, the linearized eigenvalue analysis for the sub-system of 
intensity variable (see below) shows that the solutions are sta- 
ble for |.Eo| 2 below a second threshold, \E \ 2 < \E t hr,2\ 2 , 
where 



\E, 



thr,2\ 



lh 2 



_7 

7o 



\E, 



thr,l\ 



J_ 

7o, 



(3.12) 
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This implies that the first above threshold regime is restricted 
to: 



\E, 



thr,l 



< \Eof < \E> 



thr,2\ 



(3.13) 



This is shown in Fig. |2] along with the steady state solutions 
ofEq. dlllO}. 

In this regime, the first nonlinear crystal operates in the 
above-threshold (stimulated) regime, while the operation of 
the second nonlinear crystal is in the below-threshold (spon- 
taneous) regime. The steady state solutions for the uj , lo\ and 
c^2 modes are the same as in nondegenerate parametric down- 
conversion with a single crystal 13011 . except that the stability 
region has now an upper bound. 



3. Second above-threshold regime 

In the second above-threshold regime, both nonlinear crys- 
tals operate with stimulated emission, and the amplitudes of 
all intracavity modes are nonzero. The mode w 2 acts as the 
pump mode with respect to the second nonlinear crystal and 
its intensity is above the respective threshold for stimulated 
down-conversion ui2 —> 0J3 + UJ4. Note that \E t hr,2\ is very 
close to |-Ethr,i| in the case of a strongly damped or nonres- 
onant primary pump ao that we consider here. This makes 
this second above-threshold regime quite accessible experi- 
mentally and, in the limit 72 — > 0, could bring about effective 
%( 3 ) behavior (see next section). 

Again using the intensity and phase variables, the steady 
state solutions can be written as follows: 



\Ea 



1 — 2 " 

(70 + 7) 



In addition, we show in the next section that the linearized 
eigenvalue analysis reveals that the sub-system of intensity 
variables is stable in this region. Thus, the second above- 
threshold regime corresponds to Eq. ( 13.17b and is pictured in 
Fig.|2]with its corresponding steady state solutions. 



IV. STABILITY PROPERTIES 

Here we give the details of the linearized eigenvalue analy- 
sis to determine stability of the classical steady-state regimes. 
In order to explain this approach, we proceed with a dimen- 
sionless analysis, in terms of a small parameter 



9 



X 
1 



(4.1) 



We now wish to derive the leading order behavior of the 
stochastic fluctuations in each mode, as an expansion in terms 
of g. It is simplest to first transform to dimensionless parame- 
ters, defining dimensionless time as: 



T = "ft 



(4.2) 



This scaled time variable will be used for all derivatives de- 
fined in this section. Furthermore we will also use the dimen- 
sionless parameter 



_ To 
7r 7 



(4.3) 



We note here that a linearized analysis is only the first stage 
in a stochastic diagram perturbation expansion! 31], which 
in general needs to be taken to higher order to reveal non- 
Gaussian behavior l 1 311 . The details of this will be treated else- 
where. 



\Eo\' 



4 (70 + if 



(3.14) 



A. Positive-P method 



4- 



4 = 0. 



(3.15) 



The intensities n\, n®, and n§ (rig = u°) are related by a 
simple relationship 

n a l =n% + nl. (3.16) 

This reflects the photon number conservation in the second 
crystal and the correlation between the photons uj\ and lo 2 , 
including the possibility of conversion of photons U2 into a 
pair of photons UJ3 and u>±. From the expressions for n® and 
n% we see that physical solutions are realized for driving field 
intensities above the second threshold, 



l-Efll 2 > \Ethr,2\ ■ 



(3.17) 



We start by using the full positive-P method to treat this sys- 
tem, together with an appropriate scaling for the below thresh- 
old fields, by introducing: 



0o 
03 
04 



(a - a°) 1 9 



Using the semi-classical steady state solutions, Eq. d2.5t , and 
dropping higher-order terms of order y/g or higher, we get: 



A) = 

&i = 

a 2 = 

k = 

Pi = 



~7r/3o - ai«2 
— a\ + ea? + ^t]i(t) 
-a 2 + ta\ + s/er) 2 {T) 
-03 + Va^V3( T ) 
-04 + Va2m( T ) 
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together with the Hermitian conjugate equations. The nonzero 
steady-state correlations of the noise terms are: 



(Vi(r)V2(r')) = S(r-r'), 
(%(t)?74(t')) = 6(t - T ), 



(4.4) 



The linearized equations for /?o and ^3 4 are all decoupled and 
have negative eigenvalues —70 and —7, respectively. Accord- 
ingly, the corresponding steady states are stable. The deter- 
ministic part of linearized equations for the remaining vari- 
ables, a\ and 02 (together with a\ and a 2 ), can be written 
in the matrix form as follows: 



a 



A u a, 



(4.5) 



where a 
given by 



(ai, a 2 , a.i , a 2 ) T and the drift matrix A is 



A" 




(4.6) 



The eigenvalues of the matrix A can be calculated explicitly, 
with the result that their real parts are all negative if 



£ 2 <1. 



(4.7) 



This defines the stability region, Eq. ( 13.8b , for the steady states 
(l3~6l l and the first threshold, Eq. ( f3T9b . 



B. Above-threshold stability 



and treated exactly without the use of linearization. This can 
be achieved by means of transforming to a new set of stochas- 
tic variables. In doing so, we note that the stochastic equations 
of motion for this system, Eqs. (12.5b . are equivalent in either 
ltd or Stratonovich formulation of the stochastic calculus. We 
employ the Stratonovich formulation which has the advantage 
that the variable changes are achieved using the usual calculus 
rules, without any extra variable-change terms. Accordingly, 
we first transform to new intensity and phase variables for the 
modes cuq, uj\, and uj 2 : 



Hj = o,o ; 



lln [£* 

2i \ 0+ 



(.7 = 0,1,2), 



(4.8) 



which we note are complex. The stochastic variables (23,4 and 

med t 

02/2 



a^" 4 , on the other hand, are transformed to 



«3 : 4 = a 3i 4e 

"3.4 ~ a 3,4 e 



(4.9) 



In these new variables, the stochastic differential equations 
become: 

h = -27 n + 2 \E \ cos (ip - O ) 

- 2xV"o«i«2 cos (0 O - 0+) , (4.10) 

hi = -27711 + 2x\M)«-i"2 cos (0 O - <j> + ) 

+ F l {t), (4.11) 

h 2 = - 27712 + 2xv / n nin 2 cos (0 O - <j> + ) 

-xV^(^i + a+a+) + F 2 (t), (4.12) 



In this section we analyze the stability of the above- 
threshold regimes. For reasons of length, we do not give a 
complete analysis of the fluctuations, but rather we simply 
determine which are the stable regimes. This allows us to 
build a complete large-signal phase-diagram, which is highly 
useful for determining the down-conversion properties of the 
cascaded device. Detailed spectral properties will be analyzed 
elsewhere. 



1. First above-threshold regime 

Inspecting the semi-classical steady state solutions, Eqs. 
(13.10b and (13. lib , we immediately notice that while the sum of 
the steady state phases 0° + 0° of the signal and idler modes is 
well defined and is equal to the phase of the driving field, ip , 
the individual values of 05 and (f> 2 remain unknown. In other 
words, there is no unique solution for the individual phases 
0" and (fy and any attempt to perform linearization around 
any chosen value of 4>\ or 0§ will generate a zero eigenvalue, 
implying that the steady states are unstable. This problem is 
known as phase diffusion l30l l32ll . 

In order to correctly analyze the set of coupled equations of 
motion in this regime, it is helpful to factorize them into a sub- 
set that can be linearized and is stable, while the equation as- 
sociated with the zero eigenvalue must be isolated (decoupled) 



; \Eo\ . 

6 = sin (0o 



nim^ (0 O _0 + ). 
no 



' 1 = Y'\/^^sin(0 o - 0+) +/i(t), 

7ll 



n ni . 
<>2 = X\l sin (0o 



"2 



X 



2iJn 2 



(0:30:4 — o 3 o 4 



o 3 = -7« 3 +XVn 2 aJ - — -J a 3 sm(0 o - 0+) 

2 V n 2 



+ . -7= (a3«4 - at a t ) a 3 + Fait), 



o 4 = -704 + XV^aJ - — \ a 4 sin (0 O - 0+ ) 

2 y n 2 



+ 



X 



4V^2 



(0304 — a|a|) 04 + Ti{t). 



(4.13) 
(4.14) 

(4.15) 



(4.16) 



(4.17) 



together with the equations for at, and oj. Here, we have 
defined the sum of the phase variables 0i and 02 via 

+ =0!+0 2 , (4.18) 
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and we note that the equations for 013 and ct4 contain terms 
that come from the time derivative of <j) 2 which have been sub- 
stituted with the right-hand side of Eq. ( 14. 151 ). 

The new noise terms in the above set of equations of motion 
are defined according to: 



F i,2 = a£ 2 \/X«oCi,2 + aiayxatCt^ 



fx, 



F3A = Vxo2~e 



Ci 



x<4 



2ia X,2 
^3,4 — 



Cl,2> 



(4.19) 
(4.20) 
(4.21) 



These must be rewritten in terms of the intensity and phase 
variables rij and <fii, for self-consistency: 



F1.2 = ^xnx,2( n o) 1/4 



-i4>i,2+i<fro/2 



(4.22) 



fx,2 = 



,1/4 



,— i<pi,2+iipo/2 



,i<t>l,2—i4>o/2^+ 



Ci 



2 ' 



^3,4 = VX("2) 1/4 C 3 ,4 



1013,4 



h 



(4.23) 



(4.24) 



We next introduce the phase sum and difference variables, 

0±=0i±0 2 , (4.25) 
and convert the equations of motion for 4>x an d 4>2 into: 



' : \v "" u/^ + v/~ ) sin (, '" ) - " j < 1 



X 



(a 3 a 4 -a + a+) + /+(<), (4.26) 

+) 



/n 2 /ni 

X 



+ W5 3 5 4 - 5+5+) + /-(*)> (4.27) 



(4.28) 



where the noise terms are 

f± = fx± /a. 



We now immediately see, that the equations of motion for 
the variables n , ni, 712, 53.4, S34, 4>o and + are decou- 
pled from the equation of motion for the phase-difference vari- 
able <f>_. All these variables except 0_ have a unique semi- 
classical steady state solution given by Eqs. (13.1 Qb - (I3 .lib , 
with (fp + = (f>i+(f>2 = (po and 5° 4 = (along with 
(5° 4 )+ = (5° 4 )* = 0). As we will show below, the lin- 
earized equations for this subsystem of variables are stable, 



and therefore these variables can be treated by means of lin- 
earization around their semi-classical steady states. Indeed, 
by introducing small fluctuations around the steady states 



Sn 0> 1.2(0 


= ^0,1,2(0 


- 77° 
"0,1, 2> 


(4.29) 


Sa 3i4 ,(t) 


= 23,4(4) - 


"3,4; 


(4.30) 


Sa+ A (t) 


= «3>(*) - 


(S3U)*, 


(4.31) 


tyo,+ (t) 


= &,+(*) - 


<+> 


(4.32) 



we obtain the following set of linearized equations: 

5n = -7o<5r7 - ~fSn + , (4.33) 

(4.34) 



5n + = ^5n + Fl(t), 



7 



-7o<50o 



2 n 
X n\ | 

7 



ty + = -2xty + + 2 7 <tyo + /+(<), 

(55 3 = -7<5a 3 + x\f^\^t + •^3 I (*)) 
554 = —7^54 + n°(55 + -I- JF4 (i), 



(4.35) 

(4.36) 
(4.37) 

(4.38) 
(4.39) 



together with the equations for <5a + 4 . Here, we have used 
the explicit expression for the steady state solution from 
Eq. (13.10b and the fact that n 4 = n°. The nonzero steady-state 
correlations of the noise terms, in the small-noise approxima- 
tion, are given by: 

(Fl(t)Fl(t')) = -(F _(t)F°_(t')) 

= 4r/n^S{t-t'), (4.40) 



(fl(t)f° + (t'))=-^6(t-t'), 



(4.41) 
(4.42) 



By substituting the steady state intensity n\ from Eqs. <I3 . 1 0b . 
the linearized equations and hence their solutions can be ex- 
pressed in terms of the driving field intensity |i?o| 2 - 

The eigenvalue analysis of the deterministic drift terms 
of the linearized equations reveals that the equations for 
Sno,+, <5n_, and (50o.+ are stable everywhere (the eigenval- 
ues have negative real parts), while the subsystem of variables 
(<55 3 , (554, (55;j~ , <55 + ) is stable only if 



\E f < # 

X 



1 + 1 

7o 



(4.43) 



This defines the second threshold, Eq. d3 . 1 2b , and hence the 
upper bound on the driving field intensity |i?o| 2 for the first 
above-threshold region, Eq. ( 13.131 ). 

The remaining equation for the phase difference variable 
4>-, Eq. ( 14.271 ), can not be linearized since the steady state so- 
lution is not well defined and linearization around any chosen 
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value 0°_ will reveal a zero eigenvalue, implying that the equa- 
tion is not stable. The right hand side of Eq. ( 14.271 ) can, how- 
ever, be simplified since all variables here can be linearized 
around their stable steady states. Thus, expanding these in 
terms of the stable steady states plus small fluctuations and 
keeping only the linear terms, we see that the deterministic 
terms all cancel each other. The resulting equation is 



0-=/°(i), 



(4.44) 



with the following nonzero correlation of the noise term: 



(/!(t)/!(f)) = 3^-t')- 



(4.45) 



equations for the intensities: 



2 X \E \ 2 X 2 

ni = —Z^ni H \Jn\n 2 cos tii ni?i2 

7o 7o 



Fi(t), 



(4.48) 



2 X \Ep\ 2 X 2 

n 2 = -2-yn 2 H V n i n 2 cosfli nin 2 

7o 7o 



2xV n 2«3^4 cos 02 + -Fa (*) > 



(4.49) 



n 3 = -2 7 n 3 + 2xV«2^3^4 cos0 2 + F 3 (t), (4.50) 



Thus, we have isolated the instability associated with a zero 
eigenvalue into a single phase variable, which is the phase 
difference 0_ between the signal and idler phases. Unlike the 
other variables, the phase difference 0_ is not a small fluctu- 
ation around a stable steady state. Instead it undergoes con- 
tinuous phase diffusion, governed by the noise term (t) in 
Eq. (14441 . 

Despite the fact that the noise terms F\ t 2 an d /i,2 ( an d 
hence and /+,-) depend explicitly on the individual 

phases of the signal and idler modes, (which are not well- 
defined), nevertheless, upon calculating the steady state noise 
correlations, Eqs. (14.401 i - (14.421 i. these phases combine into 
the phase sum + = 0i+ 02 which has a well defined steady 
state value and is stable. As a result, calculation of observ- 
ables via the solutions of the linearized equations of motion, 
Eqs. ( 14.331 ) - ( |4.39t which ultimately depend on the noise cor- 
relations - is a well defined procedure, and is independent on 
the individual phases 0i and 2 . 



C. Second above-threshold regime 

In the second above threshold regime, both parametric 
down-converters operate in the above-threshold regime. In 
addition to the phase diffusion in the signal and idler modes 
ujx and uj 2 , we now have a second source of instability which 
comes from the phase diffusion in the secondary signal-idler 
modes, lo 3 and u>±. To simplify our analysis, we assume here 
that the damping constant of the pump mode 70 is much larger 
than the damping constants of all the other modes, 



70 > 7- 



(4.46) 



Under this condition, one can adiabatically eliminate the 
pump mode from the equations of motion, Eq. d2.5l ). and re- 
strict ourselves to the dynamics of the remaining modes w\, 
L02, W3, and U14. Thus, we assume that do = during the 
evolution of the amplitudes 1x1,2,3,4, and we use the resulting 
expression for a , 



a 



1 

7o 



(E - xai"2) , 



(4.47) 



(together with the expression for a[j~) in the equations for 
a i,2,3 4- Transforming then to the intensity and phase vari- 
ables, as in Eq. (14.81 l. we obtain the following set of stochastic 



fi4 = -2jTi4 + 2xy / n 2 n 3 n4 cos ti 2 + Fi(t), (4.51) 
Here, we have defined 



h = 01 + 02 - fo, 

h = 03 + 04 - 02 ■ 



(4.52) 
(4.53) 



which can serve as a new pair of phase variables, traded in 
favor of the the signal and idler phases <f>i and 2 - 

The stochastic equations of motion for the phase variables, 
which we write at once in terms of Ox, 9 2 , 03 and 04, are: 



x\Eo\ 

7o 



sin 6*1 



X\ sin 2 + f 01 (* , 

n 2 



(4.54) 



n 2 n 3 



+ 



714 

x\E \ [m n . 
70 



n 2 n4 
n 3 



n 3 U4 
n 2 



sin t) 2 



— sin 0i + fe 2 (t), 
n 2 



n 2 n 4 . 
-X\l sin 2 + f 3 (t), 

n 3 



n 2 n 3 

-X\ S11102 +/ 4 t). 

714 



In the above equations, the noise terms are given by 

f$i = fi + h, 
fe 2 = h + fi — /2, 



(4.55) 
(4.56) 
(4.57) 



(4.58) 
(4.59) 



and 



' — {Eq - X«i"2)Ci,2 



2iai )2 V 70 



2iaJ V 7o 



_ VX«2 
73,4 — 7T- S3,4 

2ia 3 4 



2ia 



x<4 

+ ^3,4- 



(4.60) 



(4.61) 



3,4 



In addition, the noise terms -^1,2,3,4 in Eqs. ( I4.48I )-( I4.51I ) are 
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given by 



F 12 (t) = a{ 2 ,/ — (E - xaiQ!2)Ci,2 
7o 



7o 



F 3 ,4(t) = aj A y/x^C 3 A + &3A\/ Xa-t (3,4- 



(4.62) 
(4.63) 



In all these noise terms the amplitude variables have to be 
expressed in terms of the intensity and phase variables for self- 
consistency. 

By inspecting Eqs. (|4.48I >-( |43T| > and Eqs. ( 14.54l >-( |4"37j >, we 
see that the equations for the intensities ^12,3,4 and phases 
#1,2 are decoupled from the equations for the phase variables, 
03,4. The variables ^12, 3, 4 and 6*1,2 all have well defined 
semi-classical steady states, c.f., Eqs. ( 13.14b and (13.15b , with 
^1 2 — 0' an d as we w iH show below, the linearized eigen- 
value analysis indicates their stability. Thus, this subsystem 
of variables can be treated within the linearized treatment of 
fluctuations. The phase variables 03 and 04, on the other hand, 
do not have stable semi-classical steady states and can not be 
treated by means of linearization. To demonstrate the stability 
of the intensities ni.2,3.4 and phases 6*1.2, we introduce fluctu- 
ations around the semi-classical steady states, 



5n j {t)=n j {t)-n% ( j = 1,2,3,4), 



,(<)-<* 





l,2i 



(4.64) 
(4.65) 



and derive the following linearized equations for the intensity 
fluctuations: 



7 



Srii = —7 1 H Sni 

7o, 



X 2 "i 

7 



1_ Z ) 5n 2 +F?(t), (4.66) 
70/ 



Sn 2 =-7(1+ — ) Sn 2 

707/ 



- j I -Z)s ni -^8n + + F°{t). (4.67) 



7 

5n_ = -2 7 (5n_ + F°_(t), 



(4.68) 
(4.69) 



where we have additionally transformed to the intensity sum 
and difference variables 



5n + = (5tt,3 + (5n.4, 
(5n_ = 5ri3 — Sn^, 



(4.70) 



to further simplify the eigenvalue analysis. We have also de- 
fined 



F^ = F°±F°. 



(4.71) 



As we see, the equation for the intensity difference <5n_ fluc- 
tuation is decoupled and immediately results in a negative 



eigenvalue in the drift term, implying stability. The coupled 
equations for fluctuations in Sni, 5n 2 , and 5n + result in a 
cubic equation for the eigenvalues of the respective drift ma- 
trix. While this cannot be solved explicitly, however, the neg- 
ative real parts of the eigenvalues required for stability is as- 
certained here using the Routh-Hurwitz criterion lf33ll . 

The linearized equations for the phase fluctuations (56*1,2 (i) 
are: 



56 l = - 



x\Eo\ 



To 








/ ng + V 










) 56 2 







(4.72) 



$89! + fl 



0-2- 



(4.73) 



The eigenvalues of the corresponding drift matrix can be 
found explicitly, with the result that they all have negative 
real parts and therefore the equations are stable. The nonzero 
steady state correlations of the noise terms in the linearized 
Eqs. (l4~66l - d4~69T > and Eqs. d4/72l -( l4~73l l are: 



(F?(t)F°(t')) = 2 7 n° 1 6(t-t'), 
(F°(t)Fl(t>)) = -(F_l(t)F _(t>)) 
= -2(F°(t)F°(t>)) 
= 4-yn%6[t-t'), 
(fl(t)fl(t')) = -2(fl(t)fl(t')) 
= 2</°(t)/°(f)) 

0. 



(flit) f 1(f)) 



7 

2</3°W/4°(i')) 

-^S(t-t>). 
nl 



(4.74) 



(4.75) 



(4.76) 



(4.77) 



The remaining phase variables, 03 and 04, cannot be 
treated within the linearized fluctuation treatment, however, 
the right hand sides of the corresponding equations of motion, 
Eqs. ( 14.561 ) and ( 14.571 ). can be simplified since all variables 
here have stable steady states and can be linearized. This gives 



63 = -7(56*2 

64 = -7(56*2 



/ 3 °(*), 
flit)- 



(4.78) 
(4.79) 



where the nonzero steady state correlation of the noise terms 
is given in Eq. ( 14.77b . To further simplify the analysis we 
introduce the sum and difference phase variables 



6*± = 03 ± 04, 

for which the equations of motions are 

e + = -2jse 2 + fl(t), 
0- = /£_(*), 



(4.80) 

(4.81) 
(4.82) 
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The source of instability for the phase variable 6L is obvious, 
while for the phase variable 9 + the presence of a zero eigen- 
value is revealed when the corresponding linearized equation 
is combined with Eqs. (l4-.72t - d4.73b . Thus the variables 6 + 
and 0_ can not be linearized, and have to be treated exactly. 
The nonzero correlations of the noise terms ffj = f$ + 
and f§_ = /° - /° are 



fe + (t)fe + (t')) = -(feAt)feAt') 
= 2(/ 3 (i)/°(i')) 



"3 



(4.83) 



From Eq. ( 14.811 1 we see that the dynamics of the phase vari- 
able 9 + depends on phase fluctuations in 662, and therefore 
the equation for 9 + has to be integrated after solving for 89 2 , 
Eqs. (14.72l )- (14.73l ). The solution for 9 + {t) can be written as 



9+{t) = 9+(t ) + / -27*02(0 + /<9 + (0 < < 4 ' 84 ) 

Jt 1 J 

while the solution for 9- (t) is 

M*) = M*o) + [ feAt'W- (4.85) 



Since <5#2(i) as a solution to the set of linearized Eqs. (14.72b - 
(|4.731 > depends on the noise terms f$ and fg, the calculation 
of correlations involving the phase sum variable 9 + (t) will 
also depend on the following nonzero noise correlation: 

(fl(t)fl(t')) =2(f°(t)f 4 ( t i)) = -^S(t-t'), (4.86) 



while {flWl{t>)) = 0. 

This completes the analysis of the system in the second 
above-threshold regime. 

V. NUMERICAL SIMULATIONS 



are integrated numerically using a fourth order Runge-Kutta 
routine similar to the method given in II22I1 . The first and sec- 
ond nonlinearities were taken to be equal, i.e., \i — X2- All 
down-converted fields were given minute initial amplitudes 
and random initial phases, of which the subsequent dynami- 
cal phases were independent. 
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Figure 3: OPO intracavity powers for 72 = 0.08,71,3,4 
0.14,7o = 2.0, \Eo\/\E t hresh,2\ = 3.5, for zero detunings. 
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A qualitative reasoning identifies the far-above-second- 
threshold situation as interesting for mimicking a x' 3 ' OPO, 
in the regime where losses for the intermediate pump S 2 are 
negligible, i.e., 72 <C 71 = 73 = 74 <C 70 (note that this 
is different to the condition given in Eq. d3.41 >). Indeed, these 
hypotheses should yield close-to-ideal down-conversion rate 
from field S2 to signal fields S3 and S4, comparable to the 
emission rate into Si, and therefore be consistent with the ex- 
pectation of threefold quantum correlations between Si, S3, 
CL4, which should be non-Gaussian (another favorable situa- 
tion for this effect would be the case \2 3> Xi)- The goal of 
the following numerical simulations is therefore to ascertain 
the stability of the resonant cascade in such cases, which are 
not covered by the previous analytical treatment. 

The numerical treatment is limited to the classical equations 
of motion, given by Eq. ( 12.5b with Q — 0, i E [0, 4], which 



Figure 4: Individual field phases, going through the second thresh- 
old, for 72 = 0.08,71,3,4 = 0.14,73 = 1.0,\E \/\E thresh>2 \ = 
3.5, for zero detunings. 



A. Steady-state solutions 

As mentioned above, we restrict our analysis to a set of 
parameters such that 72 < 71 = 73 = 74 *C 70- That 
is, the primary pump mode is not resonant, the intermediate 
pump is highly resonant (most of its losses occur in down- 
conversion), and the signal fields are sufficiently resonant to 
acquire a threshold as low as a typical single-stage, doubly 
resonant OPO (DRO). In this case the OPO fields show de- 
caying oscillations to a steady state after reaching the second 
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Figure 5: The phase differences 9i and 02 between the fields in- 
volved in the first and second stages of the nondegenerate interac- 
tion, going through the second threshold, for 72 = 0.08,71,3,4 = 
0.14,73 = 1.0, \E \/\E t hresh,2\ = 3.5, for zero detunings. The 
stationary phase differences above the second threshold indicate that 
down-conversion is taking place in both stages. 
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Figure 7: Top: OPO intracavity powers for (70,72,71) = 
(10,0,0.02), \E Q \/\E thresh , 2 \ = 0.7. Bottom: Phase differences 
0i and 82, for the same conditions (phases are numerically wrapped 
inside [— tt, it]). 



B. Second above-threshold regime for low 72 



Figure 6: OPO intracavity powers for 72 = 0.001,71,3,4 = 
0.14,7o = 1.0, \Eo\/\E t hresh,2\ = 5.8, for zero detunings. Note 
that «4 is omitted since its plot follows exactly 0:3. The spiking fre- 
quency increases with l-Eol/l-Ethresh^l when other parameters are 
held constant. 



threshold. Higher losses for 72 result in over-damping. Fig- 
ure[3]illustrates the steady-state solutions. Both thresholds are 
clearly visible. At long times (t > 100 in Fig. |3) the field 
amplitudes match the stationary solutions of the previous sec- 
tion. Increasing Eq or 72 causes the oscillations to decay more 
quickly. 

Figure |4] shows the individual phases as the second thresh- 
old is reached. The final phase in the steady state is inde- 
pendent of the initial starting phase of any of the fields' seed 
values. (The primary pump parameter Eq is taken to be real.) 
Figure [5] shows the nonlinear phase differences 9\ and #2 for 
the first and second stages of the OPO. It is clear from these 
phase differences that the system is in a state of cascaded para- 
metric down-conversion. 



1. Nondegenerate cascade 

One interesting question is that of obtaining a stable effec- 
tive x OPO by lowering 72 and operating well above sec- 
ond threshold. In that case, the nondegenerate and degenerate 
cascades do not exhibit the same behavior. As the intermedi- 
ate pump loss rate is lowered, a spiking instability is obtained 
in both cases, as displayed in Fig. [6] for the nondegenerate 
case. One may overcome this self-pulsing and induce a tran- 
sition to a stable steady state by increasing the pump parame- 
ter above threshold. The lower 72 the higher l-Eol/l-EWesMl 
needs to be to achieve steady state in the nondegenerate case. 



2. Degenerate cascade 

In the degenerate case, a remarkable result is that the spik- 
ing behavior is always transient and relaxes into a stationary 
state. However, more insight into the behavior of the degener- 
ate cascade is obtained, once again, by scrutinizing the evolu- 
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Figure 8: Phase differences 8i and 62, for zero detunings and 
\E a \/\E th resh,2\ = 2.0. Top: 72 = 0.08, 71 = 0.14. Bottom: 
71 = 72 = 0.14 (phases are numerically wrapped inside [—it, n]). 



Figure 9: OPO intracavity powers for 70 = 
0.14, \E \/\E t hresh,2\ = 2.8, for zero detunings. 
Bottom: 72 = 0.14. 



0.08,71,3,4 = 
Top: 72 = 0.09. 



tion of the phases Q\ and 62- This bears particular physical sig- 
nificance for the degenerate cascade because the first stage can 
stimulate emission in the second stage, which cannot happen 
in the nondegenerate cascade due to the indistinguishability 
of the signal fields 13411 . Because of this effect, the degenerate 
cascade will exhibit greater sensitivity to the evolution of d\ 
and 6*2, whose swings translate into the appearance of compet- 
ing sum-frequency generation (SFG) processes in both stages. 

a. First above-threshold: Figure [JJ (top) displays the 
damping of the low 72 spiking. The phase behavior is plotted 
in Fig. [JJ (bottom) and shows that significantly rich transient 
evolution eventually yields a true cascade of two parametric 
down-conversions (PDCs). 

b. Second above-threshold: In Fig. [8] (top), PDC is not 
the only process occurring in both stages: the solutions can be 
seen to have a PDC component and a competing SFG com- 
ponent. This is consistent with the entering of the stimulated 
emission regime in the second stage as one crosses the sec- 
ond threshold. The system is able to find a steady state so- 
lution nonetheless. However, the quantum statistics might be 
expected to be nontrivially affected. If one increases 72 to 
the level of the signal loss rate, the phase evolution yields 
this time a stable PDC cascade (Fig. [8] (bottom)), as already 
demonstrated analytically in the previous section. 



In conclusion, the degenerate cascade, because of the ad- 
ditional signal feedback between the two stages, is clearly a 
much richer system than the nondegenerate cascade. This ad- 
ditional feedback leads to stabilization of the PDC cascade in 
the low 72 regime in the first above-threshold regime. In the 
second above-threshold regime, however the degenerate cas- 
cade displays two stable regimes, one of which does not have 
pure PDC character. Bistable behavior or a bifurcation is to 
be expected there. This also opens interesting horizons for the 
quantum simulations of degenerate resonant cascades. 

C. Second above-threshold regime for low 72 , 70 

It is interesting to briefly investigate the behaviour found 
for low 70. This regime involves a low-loss, resonant pump 
mode. It has different stability properties to the situations 
treated elsewhere in this paper, and more dramatic behav- 
ior is observed. By setting 70 ~ 72 < 71,3,4, c.f. Fig. [9] 
(top), the OPO becomes unstable above the second thresh- 
old, where the pump ciq and intermediate pump a-2 compete 
strongly. The above-threshold phase leads to a return to the 
first-above threshold regime, before recurring. In the case 
where 70 < 72, c.f., Fig. [9] (bottom), the amplitude of the 
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Figure 10: Evolution of intracavity fields for 70 = 1.0,72 = 
0.08,71,3,4 = 0.14, \E \/\E threah>2 \ = 3.5, perturbation S ai = 04, 
for zero detunings. The perturbation has been applied simultaneously 
to the real and imaginary parts of the fields, leaving no net perturba- 
tion to the phases. 



Figure 12: Evolution of intracavity fields for 70 = 1.0,72 = 
0.08, 71,3,4 = 0.14, \Eo\/\E threah>2 \ = 3.5, perturbation S ai = a„ 
for zero detunings. The perturbation has been applied to the real part 
of the field amplitudes. 
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Figure 11: Evolution of the phase differences 9i and 62 for 70 = 
1.0,72 = 0.08,71,3,4 = 0.14, \Eo\/\E t hreshfl\ = 3.5, perturbation 
Sai = di, for zero detunings. The perturbation has been applied to 
the imaginary parts of the fields. 

oscillations above the second threshold keep increasing and 
the system never reverts to the first-above threshold regime. 



For small perturbations on the order of a couple percent of 
the steady state amplitudes (perturbations to phase included), 
the system returns to steady state after a few oscillations, and 
the field phases also return to their steady state values. Some 
perturbations may change the individual steady-state phases; 
however these changes are inconsequential if the compound 
phase differences 6-y and 82 remain at the PDC values. Fig- 
ure [TO] shows a typical response of the intracavity powers to 
a perturbation on the field amplitudes only. The OPO returns 
to the original steady state solutions even under quite large 
disturbances. Figure [TJJ shows an important part of the phase 
evolution. Upon disturbing the system, the phase of the sec- 
ondary pump shifts by ir and then back by — ir. This is only 
true for a large change in the real or imaginary components of 
the field (greater than 50% in this case). The phases of all the 
other fields remain comparatively unaffected. Thus, the phase 
differences 6\ and 82 shift by ir quickly and then by — ir (see 
Fig. [TlTl, When the disturbance is small ~ 10%, the phase dif- 
ferences recover so quickly that a change in phase differences 
is not observed. 



D. Stability analysis of stationary solutions 

We simulate the effect of a perturbation by causing an in- 
stantaneous change in the intracavity fields and observing the 
numerical response of the system. Of particular interest is the 
phase evolution of the stationary solutions under two different 
types of perturbations, for this will give insight into compet- 
ing interaction (SFG/PDC) behaviors in the degenerate case. 
We distinguish several types of perturbations: (i) amplitude 
changes of the fields, leaving the phase unperturbed. ( ii) phase 
change of the fields. (Hi) change in both. We display typical 
results obtained for a variety of magnitudes of change. 



Figure [12] shows a more complicated amplitude response 
when a perturbation is applied to both the phase and ampli- 
tude simultaneously. The OPO recovers the stationary ampli- 
tudes, but each field, except the primary pump, also undergoes 
a permanent phase change, c.f. Fig. [13] (top), even though 9\ 
and 82 return to unaltered values after opposite fluctuations, 
c.f. Fig. [13] (bottom). With increasing pump parameter, the 
phase changes in Fig. [13] and Fig. Q~T] can undergo several 
sign changes before returning to steady state. This effect can 
be seen when the applied perturbation is very large. These 
results indicate that nondegenerate cascade is essentially as 
stable as a single-stage DRO. 
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Figure 13: Evolution of the OPO phases for 70 = 1.0,72 = 
0.08,71,3,4 = 0.14, \E \/\E t hresh,2\ = 3.5, perturbation 5 ai = a,:, 
for zero detunings. The perturbation has been applied to the real part 
of the fields. Top: individual phases. Bottom: 9i and 62 (phases are 
numerically wrapped inside [— ir, 7r]). 

VI. CONCLUSION 

In conclusion, we have given a preliminary analysis of the 
novel properties of a doubly cascaded nondegenerate intra- 
cavity parametric oscillator. This has the property that it 
is able to mimic a x*- 3 -* down-conversion system, while still 
relying on the properties of widely available phase-matched 



%( 2 ) down-conversion crystals. Our analysis focuses on con- 
structing phase-space equations for the cascaded system, and 
demonstrating the existence of multiple thresholds and stable 
regions. 

In the case of five non-degenerate modes, we have de- 
rived phase-space equations in the full double-dimensional 
positive-P representation, as well as approximate equations 
using the semi-classical or truncated Wigner approach. We 
show the presence of three distinct classically stable regimes, 
corresponding to below-threshold operation, an intermediate 
threshold where only some of the modes are classically ex- 
cited, and a fully above threshold regime similar to x^ 3 -* down- 
conversion. 

A detailed analysis of stability of these regimes is carried 
out to show whether the relevant driving fields will result 
in stable operation. This analysis is restricted to the non- 
degenerate case, for parameter values in which all losses are 
equal except for the pump, which is assumed to be strongly 
damped. The nonlinear coefficients are also assumed to be 
equal. We find that for these parameter values each of the 
three regimes mentioned is stable, that is, small signals are 
damped back to the classical steady-state values. 

We also give dynamical simulations of the mean field 
equations, which clearly demonstrate the existence of sta- 
ble regimes, as well as unusual phase-evolution and distinct 
dynamical behaviour in the degenerate and non-degenerate 
cases. Remarkable coincidences of two [35] and even three 
lf36h nonlinear interactions in a single-grating periodically 
poled crystal have been observed, which illustrates the ex- 
perimental possibilities of such a technique. The dynamical 
analysis in this case, although based on classical equations, 
is able to treat a larger variety of parameters and detunings, 
as well as allowing an investigation of stability in the case of 
much larger perturbations. The general conclusion is that both 
the cascaded DPO and NDPO have a rich variety of stable op- 
erating regimes and thresholds, including the possibility of an 
above threshold x^ 3 -* domain. 
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